Detecting Seasonality in individual indicators
load(paste0(IO$output_clue,"indiv_indicators.Rdata"), verbose = TRUE)
## Loading objects:
## indiv_indicators
colnames(indiv_indicators)
## [1] "cycle_id" "user_id"
## [3] "date" "country"
## [5] "age_cat" "bmi_cat"
## [7] "active_use" "day_id"
## [9] "bleeding" "cycle_nb"
## [11] "cycle_length" "cycleday"
## [13] "cycleday_from_end" "ovu_day"
## [15] "daily_fertility" "time_num"
## [17] "mucus" "mucus_num"
## [19] "sd_cycle_length" "n_eggwhite"
## [21] "mucus_quality" "n_bleeding_after_day_5"
## [23] "n_heavy_bleeding" "cycle_fertility"
## [25] "age" "bmi"
## [27] "overall_fertility" "fertility"
for(user in unique(indiv_indicators$user_id)[1:20]){
this_user_indiv_indicators = indiv_indicators[indiv_indicators$user_id == user,]
this_user_indiv_indicators = this_user_indiv_indicators[order(this_user_indiv_indicators$time_num),]
if(!all(this_user_indiv_indicators$fertility == 0)){
timeserie_matrix = as.matrix(this_user_indiv_indicators[,c("time_num","daily_fertility")])
current.wt = wt(timeserie_matrix,
dt = par$dt, dj = 15*par$dt,
lag1 = 0.7, # pink noise
sig.test = 1) # we want to run a significance test
plot(current.wt, main = paste0(user," - daily_fertility"))
timeserie_matrix = as.matrix(this_user_indiv_indicators[,c("time_num","fertility")])
current.wt = wt(timeserie_matrix,
dt = par$dt, dj = 15*par$dt,
lag1 = 0.7, # pink noise
sig.test = 1) # we want to run a significance test
plot(current.wt, main = paste0(user," - fertility"))
}
}
Detecting Seasonality in population indicators
load(paste0(IO$output_clue,"pop_indicators.Rdata"), verbose = TRUE)
## Loading objects:
## pop_indicators
agg_col_combinations = unique(pop_indicators[,par$agg_col])
for(i in 1:ncol(agg_col_combinations)){
agg_col_combinations[,i] = as.character(agg_col_combinations[,i])
}
for(i in 1:nrow(agg_col_combinations)){
cond = rep(nrow(pop_indicators), TRUE)
for(agg_col in par$agg_col){
eval(parse(text = paste0("cond = cond & (pop_indicators$",agg_col," == agg_col_combinations$",agg_col,"[i])")))
}
i_cond = which(cond)
if(length(i_cond)>0){
this_pop_indicators = pop_indicators[i_cond,]
# daily fertility
timeserie_matrix = as.matrix(this_pop_indicators[,c("time_num","daily_fertility")])
timeserie_matrix = timeserie_matrix[!is.na(timeserie_matrix[,2]),]
LAG1 = arima(timeserie_matrix[,c('daily_fertility')], order = c(1,0,0))$coef[1]
current.wt = wt(timeserie_matrix,
dt = par$dt, dj = 15*par$dt,
lag1 = LAG1, # autocorrelation
sig.test = 1) # we want to run a significance test
plot(current.wt, main = paste0("daily fertility", paste0(agg_col_combinations[i,], collapse = " | ")))
# female fertility
timeserie_matrix = as.matrix(this_pop_indicators[,c("time_num","fertility")])
timeserie_matrix = timeserie_matrix[!is.na(timeserie_matrix[,2]),]
LAG1 = arima(timeserie_matrix[,c('fertility')], order = c(1,0,0))$coef[1]
current.wt = wt(timeserie_matrix,
dt = par$dt, dj = 15*par$dt,
lag1 = 0.7, # pink noise
sig.test = 1) # we want to run a significance test
plot(current.wt, main = paste0("fertility", paste0(agg_col_combinations[i,], collapse = " | ")))
}
}
